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I.  Summary  of  Research  Activities 

This  grant  started  on  June  1,  2001,  and  in  the  following  we  report  on  research  activities 
from  that  date  to  November  30, 2003. 

Electrified  Jets  and  Electrohydrodynamics  of  molecular  fluids  at  the 
nanometer  scale. 

A  large  part  of  the  initial  phase  of  this  work  involved  the  development  and  testing  of  the 
molecular  dynamics  computer  code  and  associated  analysis  routines,  to  be  used  in  the 
study  of  electrified  jets  and  electrohydrodynamics  of  a  number  of  common  molecular  and 
ionic  fluids  that  are  currently  being  used  experimentally  in  these  areas  -  along  with 
fundamental  interest  in  these  systems,  these  studies  aim  at  the  application  of  electrified 
jets  to  colloid  thruster  technology.  A  major  portion  of  our  efforts  was  devoted  to  advance 
the  computer  algorithms  to  the  point  that  we  could  simulate  the  operation  of  a  colloid 
thruster  and,  as  will  be  discussed  below,  we  have  entered  this  latter  stage. 

Existing  studies  of  electrified  jets  and  electrohydrodynamic  phenomena 
specifically  involving  the  interaction  of  fluids  with  externally  applied  electric  fields 
utilize  continumn  modeling  and  there  have  been  no  prior  atomistic  simulations  in  these 
areas.  A  significant  component  of  our  progress  to  date  has  been  the  development  of 
parallel  computer  algorithms  that  can  efficiently  model  systems  as  complex  as  a  colloid 
thruster,  yet  be  flexible  enough  to  allow  modeling  of  other  interesting  phenomena 
involving  the  electrodynamics  of  dielectric  and  charged  fluids. 

Since  this  project  initiates  a  novel  molecular-scale  approach  to  the  investigation 
of  electrohydrodynamic  phenomena,  we  have  chosen  to  adopt  a  strategic  staged 
approach,  where  we  developed  our  research  in  stages  of  progressively  increasing 
complexity.  Our  work  falls  roughly  into  two  phases.  The  first  phase  encompasses  the 
development  of  computer  code  of  parallel  architecture  that  uses  well-tested  atomic  and 
molecular  interaction  potentials  and  most  importantly  the  simulation  code  must  compute 
in  a  very  efficient  maimer  the  long-range  Coulombic  interactions  between  charged 
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species.  The  second  phase  involves  building  upon  the  simulation  code  described  above  so 
that  one  can  model  the  complex  geometries  and  electric  fields  present  in  a  colloid 
thruster. 


Phase  One:  Long-Range  Interactions 

As  an  early  test  of  the  developing  computer  code  and  as  a  guide  in  directing  furttier  work, 
we  first  studied  single  droplets  of  formamide  with  varying  concentrations  of  sodium 
iodide  in  an  external  uniform  electric  field.  This  is  a  classic  problem  in  the  study  of 
electrified  fluids  that  has  been  studied  both  experimentally  and  theoretically  and  we 
expected  that  the  comparison  of  our  results  to  theoretical  predictions  would  be  very 
helpful.  In  fact,  much  of  the  important  physics  that  is  involved  in  colloid  thrusters 
operating  in  the  cone-jet  mode  is  already  present  in  the  simpler  problem  of  the  elongation 
of  a  conducting  liquid  droplet  in  a  uniform  external  electric  field.  The  fluids  that  we 
currently  focus  on  are  the  molecular  fluid  formamide  along  with  ionic  fluids,  particularly 
sodium  iodide  which  is  one  of  the  common  ionic  salts  added  to  formamide  to  increase  its 
conductivity  and  effectiveness  in  producing  electrified  jets  in  colloid  thrusters. 

In  order  to  simulate  systems  that  are  sufficiently  large  to  enable  us  to  make 
contact  with  continuum  behavior,  the  systems  that  we  study  could  eventually  contain 
hundreds  of  thousands  of  atoms  and  it  is  important  to  utilize  simple  efficient 
representations  of  the  inter-atomic  forces.  Formamide  is  a  planar  molecule  (see  Fig.l) 
whose  internal  degrees  of  freedom  and  intra-molecular  forces  are  of  lesser  significance 
for  the  phenomena  in  which  we  are  interested  here.  Consequei^y,  we  treat  the 
formamide  molecule  as  a  solid  body  using  quaternion  dynamics  (ref.  1),  implemented  via 
a  mid-step  implicit  leap-frog  algorithm  (ref  2)  that  has  been  shown  to  be  an  extremely 
stable  integration  scheme  (witii  only  a  very  small  energy  drift  occurring  for  long 
simulation  periods).  The  geometry  of  the  formamide  molecule  is  taken  from  high- 
resolution  X-ray  studies  of  formamide  crystals  (ref  3). 

In  our  simulations  we  employed  the  AMBER  force  field  parameters  (refs.  4-5)  for 
file  intermolecular  van  der  Waals  interactions  between  the  atomic  sites  located  on 
different  molecules,  and  also  for  interactions  between  the  formamide  molecule  and  the 
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ionic  species  (e.g.  sodium  iodide),  as  well  as  for  the  description  of  ion-ion  interactions. 
For  the  atoms  of  the  formamide  molecule  we  use  the  CHELP-BOW  (ref.  6)  partial 
charges;  these  have  been  shown  to  give  a  good  overall  description  of  the  electrostatic 
potential  of  the  molecules.  The  weak  inter-molecular  van  der  Waals  interactions  are 
truncated  on  a  group  basis  through  the  use  of  a  smooth  switching  function  (ref.  7)  that 
depends  on  the  distance  between  molecular  centers  of  mass  so  that  entire  groups  on  one 
molecule  interact  with  the  entire  group  of  atoms  on  another  molecule,  and  dipole 
interactions  are  not  truncated. 

Our  very  first  effort  was  to  develop  and  test  a  computer  code  of  parallel 
architecture  for  simulations  of  pure  formamide  and  then  we  proceed  to  add  the  ionic 
components.  As  the  computer  code  was  developed  we  performed  simple  studies  of  liquid 
and  solid  formamide  systems  to  test  the  code  and  to  assess  the  ability  of  the  interatomic 
potentials  and  the  model  that  we  use  to  reproduce  basic  known  properties  (such  as  solid 
and  liquid  densities  and  the  melting  point  of  formamide). 

Since  this  is  the  first  large-scale  atomistic  simulation  involving 
electrohydrodynamics  of  complex  fluids,  we  had  to  learn  from  our  initial  studies  what 
features  are  essential  in  our  computer  code  to  faithfully  model  such  systems.  In 
particular,  our  initial  code  utilized  long-range  Coulombic  interatomic  potentials  with 
finite  interaction  cutoffs  as  this  was  a  simple  first  approach  that  is  sufficient  in  many 
simulations.  However  when  we  ran  test  studies  of  pure  formamide  droplets  in  uniform 
external  electric  fields  of  varying  strengths,  we  noted  significant  cooperative  effects 
within  the  dipolar  fluid  that  emerged  as  the  cutoffs  were  increased  fi-om  values  of  ~  9  -  12 
A  (typical  values  in  many  studies  involving  atomic  partial  charges)  to  36  A.  Droplets 
(which  when  studied  with  small  interaction  cutoffs  remain  spherical)  elongate  as  the 
cutoff  is  increased  due  to  cooperative  dipole-dipole  and  dipole-field  interactions,  in 
qualitative  agreement  with  existing  continuum  theories  of  dielectric  fluids  in  electric 
fields  (ref  12).  From  these  observations  it  became  apparent  that  the  long-range 
interactions  need  to  be  treated  much  more  accurately.  The  necessity  of  this  became  even 
more  evident  when  our  computer  code  was  advanced  to  the  point  of  allowing  the  addition 
of  ionic  salts  to  the  dielectric  fluid  -  their  description  especially  requires  a  proper 
treatment  of  the  long-range  Coulomb  forces.  The  development  of  our  computer  code  to 
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address  these  issues  is  one  of  the  most  important  parts  of  our  progress  and  we  will 
describe  some  of  its  capabilities  in  more  detail,  as  well  as  discuss  some  of  the  issues  that 
we  have  begun  to  study  with  it. 

Fast  multipole  method  (FMM).  One  of  the  most  highly  developed  methods  to 
efficiently  model  long-range  forces  is  the  fast  multipole  method  (ref  8),  in  which  a 
system  is  spatially  represented  in  a  large  (roughly)  cubic  “root-cell”.  The  latter  contains  a 
hierarchy,  or  “tree-structure”  division,  of  “child”  sub-cells  obtained  by  dividing  the  root¬ 
cell  into  octants  and  then  allowing  these  child  cells  to  be  the  parent  cells  of  further  binary 
subdivisions  until  one  reaches  the  smallest  “leaf’  cell  (see  Fig.2).  Multipole  moments  are 
computed  for  all  cells  at  all  levels  from  the  atoms  contained  within  the  cells.  Atoms  in  the 
basic  leaf  cell  interact  directly  with  all  the  atoms  in  neighboring  leaf  cells,  while  they 
interact  with  the  multipole  moments  of  more  remote  cells;  the  more  remote  the  cells  are 
from  the  central  leaf  cell,  the  larger  tiiey  are  allowed  to  be.  Thus  an  atom  interacts  with  a 
small  set  of  multipoles  representing  entire,  increasingly  larger,  volumes  of  remote  space, 
rather  than  with  all  of  the  atoms  contained  within  it  -  this  is  the  main  feature  that  makes  it 
possible  to  efficiently  model  large  systems  involving  long-range  interactions. 

Initially,  we  searched  for  existing  FMM  computer  code  that  could  be  easily 
adapted  for  our  needs  as  the  FMM  algorithm  cSh  be  quite  complex.  While  the  basic  FMM 
technique  has  been  around  for  over  a  decade,  its  various  forms  and  implementations  are 
only  beginning  to  evolve  into  standardized  codes  that  can  be  simply  adapted  for  general 
studies.  We  did  not  find  an  existing  code  that  could  be  easily  modified  for  the  various 
types  of  systems  and  geometries  of  potential  interest  to  us.  Consequently,  we  further 
developed  our  computer  code  by  formulating  and  programming  our  own  parallel  FMM 
algorithms.  This  affords  us  a  lot  of  flexibility  in  our  ability  to  modify  and  tailor  our 
programs  for  different  studies.  Some  of  the  features  of  our  code  are  non-standard  and 
many  of  the  more  standard  features  are  not  often  found  coexisting  in  one  computer 
program. 

Some  of  the  main  capabilities  we  have  established  in  our  code  are: 

1)  The  ability  to  replicate  basic  root-cells  and  “stack”  them  along  a  particular  axis  to 
create  a  long  narrow  tree-structure  (requiring  an  appropriately  modified  non-standard 
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FMM  algorithm).  This  is  particularly  useful  for  efficient  simulations  of  systems  such  as 
highly  elongated  droplets  and  long  liquid  jets. 

2)  A  large  degree  of  flexibility  in  the  choice  of  boundary  conditions  (he’s).  We  can 
employ  periodic  boundary  conditions  (pbe’s)  in  one,  two  or  three  spatial  directions,  or 
perform  “cluster  simulations”  where  pbe’s  are  not  imposed  at  all.  In  simulations  using 
pbe’s  in  less  then  3  directions,  the  non-periodic  directions  can  be  reflecting  or  absorbing. 
Absorbing  bc’s  (any  atoms  that  cross  the  system  boundary  are  removed  fi-om  the 
simulation)  are  necessary  in  simulations  of  fragmenting  droplets  or  liquid  jets  where  there 
are  large  fluxes  of  atoms  crossing  the  finite  boundaries  (maintaining  these  atoms  of  the 
root-cell  would  cause  artificial  effects). 

3)  A  modified  FMM  algorithm  so  that  when  using  pbe’s  the  tree-structure  of  the  root¬ 
cell,  along  with  all  its  computed  multipoles,  can  be  replicated  any  number  of  (specified) 
times  as  neighboring  image  cells.  When  this  is  done,  atoms  in  the  central  computational 
cell  see  the  correct  multipole  structure  as  they  interact  with  neighboring  periodic  image 
cells  out  to  any  desired  distance. 

4)  The  ability  to  perform  constant  pressure  simulations.  The  spatial  dimensions  of  the 
root-cell  can  be  allowed  to  vary,  while  using  three-dimensional  pbe’s  and  the  FMM, 
according  to  a  modified  constant-pressure  algorithm  due  to  Berendsen  (ref  1).  This 
allows  us  to  prepare  bulk  samples  of  both  crystals  and  dielectric  fluids  with  solvated  ionic 
salts  at  any  pressure  or  temperature  of  interest.  These  systems  are  useful  for  studying 
intrinsic  bulk  properties  of  materials  currently  being  used  in  colloid  thruster  research.  It 
also  permits  us  to  readily  “carve  out”  droplets  of  various  selected  radii  and  charge  for 
furtiier  more  specialized  studies. 

As  a  simple  check  of  both  the  computer  code  and  the  quality  of  the  interatomic 
potential  parameters,  we  performed  a  study  that  aims  at  finding  the  melting  point  of  the 
simulated  formamide  for  comparison  with  experimental  dafei.  Many  experimental  studies 
of  liquid  formamide  are  performed  at  room  temperature,  which  is  only  about  25K  above 
the  experimental  melting  point  of  formamide.  The  potential  parameters  being  used  by  us 
were  not  developed  to  specifically  rq)roduce  the  melting  point.  Therefore,  in  anticipation 
of  our  simulations  that  involve  formamide  in  the  liquid  state  close  to  the  melting  point,  it 
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is  important  to  ascertain  the  thermodynamic  properties  of  the  simulated  formamide.  Our 
results  from  constant  pressure  simulations  and  liquid-solid  coexistence  studies  indicate 
that  that  our  description  of  formamide  yields  a  liquid  density  that  is  in  rather  good 
agreement  with  experiments  (to  within  a  few  percent)  and  a  melting  point  that  is 
approximately  10  degrees  K  above  the  experimental  value  of 275K. 

Dielectric  droplets  in  high  electric  fields 

With  the  above  in  mind,  we  can  now  discuss  the  types  of  simulations  being  initiated 
using  our  parallel  FMM  code.  First,  we  give  some  general  information  about  the  types  of 
systems  we  explore  and  the  materials  involved.  As  mentioned  earlier,  these  first 
simulations  involve  single  droplets  of  formamide  with  varying  concentrations  of  sodium 
iodide  in  an  external  uniform  electric  field.  Dielectric  and  charged  droplets  in  uniform 
electric  fields  are  well-studied  problems  involving  electrified  fluids.  Some  of  the  most 
relevant  physics  involved  in  colloid  thrusters  operating  in  the  cone-jet  mode  may  be  seen 
already  in  the  simpler  problem  of  an  elongated  conducting  liquid  droplet  in  a  uniform 
external  electric  field,  including  the  formation  of  electrified  jets  and  emission  and 
acceleration  of  small  charged  droplets  at  the  ends  of  the  elongated  parent  drops. 

Our  initial  studies  are  performed  at  T=310K  (see  ref.  9)  and  involve  droplets  with 
a  diameter  of  10  nm.  These  droplets  contain  -7150  formamide  molecules  (i.e.  6  x  7150  - 
43000  atoms  with  partial  charges).  The  droplets  with  added  ionic  salts  have  additional 
Nal  molecules  in  two  concentrations;  in  one  the  number  of  formamides  to  Nal  molecules 
is  8:1  (-29  %  Nal  by  weight,  -900  Nal  molecules,  i.e.  900  Na'^,  900  T  ions  solvated  in 
the  formamide  fluid),  and  in  the  other  the  number  of  formamides  to  Nal  molecules  is  16:1 
(-17  %  Nal  by  wt,  -450  Nal  molecules).  The  selection  of  these  specific  concentrations, 
as  well  as  the  choice  of  materials  (formamide  and  Nal),  were  guided  by  current 
experiments  on  electrified  jets  and  their  role  in  colloid  thrusters  (refs.  10-11,13),  in 
conjunction  with  the  desire  to  allow  comparisons  with  experimental  findings. 

It  is  noteworthy  that  the  droplets  simulated  here  are  large  enough  to  allow 
meaningful  comparisons  with  continuum  electrohydrodynamic  theories,  yet  they  are 
sufficiently  small  to  enable  us  to  perform  (in  a  reasonably  short  time  frame)  a  wide  range 
of  initial  computer  experiments.  These  studies  should  not  only  be  relevant  in  order  to 
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gain  a  better  understanding  of  electrified  jets  -  rather,  they  also  allow  a  rapid 
development  of  analysis  tools,  metiiodologies  and  basic  insights  that  will  be  important  in 
the  study  of  much  larger  systems. 

In  these  early  droplet  simulations  the  basic  root-cell  of  the  FMM  code  has  8  tree 
levels  with  the  smallest  leaf  cell  chosen  to  be  1.2  nm  on  each  side.  This  gives  a  root-cell 
of  *  1.2  nm  =153.6  nm  dimensions  on  each  side.  This  root-cell  is  replicated  3  times 
along  the  z-axis  (the  electric  field  direction)  to  give  a  total  tree-structure  of  153.6  nm  x 
153.6  nm  x  460.8  nm  in  the  center  of  which  we  place  our  10  nm  diameter  droplet  (see 
Fig.3).  The  very  large  system  size  allows  enough  room  for  the  parent  droplet  to  move  and 
elongate  freely  and  for  secondary  emissions  to  occur  jfrom  any  charged  clusters  leaving 
the  parent  drop.  Absorbing  bc’s  are  used  so  that  atoms  reaching  the  system  boundaries 
are  removed.  Once  die  simulation  begins,  no  form  of  temperature  control  or  control  of  the 
droplets  center-of-mass  motion  is  performed  as  we  do  not  wish  to  risk  introducing 
artificial  effects. 

Initially,  we  studied  a  pure  formamide  droplet,  where  a  10  nm  diameter 
formamide  droplet  is  placed  in  a  constant  uniform  electric  field  along  a  given  axis  (the  z- 
axis).  The  field  strengths  we  use  are  on  the  order  of  1  V/nm,  a  value  that  is  of  relevance 
to  current  research  in  these  areas  (ref  1 1);  this  value  of  the  field  is  large  enough  to 
generate  a  pronounced  elongation  of  the  formamide  droplet.  These  simulations  have 
begun  and  we  plan  to  vary  the  field  strengths  over  a  wide  range  and  compare  the  results 
using  these  nano-scale  dielectric  droplets  to  existing  predictions  of  continuum  theory  (ref 
12).  Our  initial  results  indicate  that  in  fields  of  -1.0-1. 5  V/nm  diese  droplets  will  exhibit 
length/width  aspect  ratios  on  the  order  of  -15:1,  (see  Fig.l),  roughly  an  order  of 
magnitude  smaller  than  the  continumn  predictions.  Additionally  we  have  discovered  that 
the  external  electric  field  can  induce  crystallization  of  a  droplet.  Namely  the  external  field 
modifies  the  thermodynamics  and  phase  behavior  of  a  dielectric  droplet.  While  these 
studies  are  not  our  main  focus,  they  also  serve  as  interesting  reference  systems  for  our 
main  study  as  they  have  no  dissolved  Nal  but  are  otherwise  similar  to  the  formamide/Nal 
droplets. 

In  parallel  with  the  simulations  of  pure  formamide  droplets  described  above,  we 
have  studied  formamide  droplets  with  dissolved  salts.  These  droplets  (10  nm  in  diameter) 
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had  several  concentrations  of  dissolved  Nal  and  they  were  placed  in  constant  uniform 
electric  fields  having  strengths  of  1.0- 1.5  V/nm.  Here  we  give  a  simple  overview  of  this 
work  and  initial  results.  The  behavior  of  these  formamide/Nal  droplets  is  significantly 
different  than  a  pure  formamide  droplet.  While  the  pure  droplet  elongates  and  reaches  a 
steady  state  aspect  ratio  of  -15:1  with  a  relatively  small  degree  of  evaporation  of 
formamide  molecules,  the  salted  droplet  (29  %  by  wt  Nal,  field  strength=1.5  V/nm) 
begins  to  deform  (within  -10-20  ps  after  the  application  of  the  electric  field)  prior  to  the 
emission  along  the  axial  field  direction,  of  positively  and  negatively  charged  clusters 
consisting  of  sodium  and  iodide  ions  solvated  by  formamide  molecules.  Over  the  course 
of  about  1  ns,  the  parent  droplet  continues  to  elongate  (to  at  most  an  -8:1  aspect  ratio)  but 
it  never  achieves  the  long  slender  end  points  seen  in  the  pure  droplet  due  to  the  constant 
stream  of  ionic  clusters  of  varying  size  breaking  fi'om  its  narrow  end  sections  and  being 
accelerated  by  the  electric  field  (see  Fig.4). 

In  Fig.5  we  display  (top  to  bottom)  an  early  stage  of  the  elongation  process, 
detail  of  negatively  charged  fragments  streaming  off  the  droplet  end  during  a  later  stage, 
and  the  system  closer  to  the  end  of  the  simulation.  Eventually,  what  remains  of  the  parrat 
drop  fragments  into  two  sizable  clusters  that  in  turn,  as  they  are  accelerated  by  the 
electric  field,  diminish  in  size  as  they  emit  a  sizeable  number  of  charged  secondary 
clusters.  When  this  salted  droplet  is  placed  in  a  lower  field  of  1.0  V/nm,  the  charged 
clusters  that  detach  are  larger.  The  smaller  field  elongates  the  droplet  more  slowly,  giving 
more  time  for  capillary  forces  to  pinch  off  (the  Raylei^  instability)  larger  sections  of  the 
droplet’s  slender  end  regions.  We  are  beginning  to  observe  and  further  explore 
phenomena  that  have  not  previously  been  studied  at  the  atomistic  level.  We  anticipated 
that  the  mechanisms  of  charged  particle  emissions  from  the  slender  ends  of  a  droplet  such 
as  the  one  studied  here,  would  be  very  similar  in  nature  to  those  emerging  from  the 
extremely  slender  jet  at  the  apex  of  a  cone-jet  in  a  colloid  thruster  -  as  we  discuss  below, 
we  have  indeed  found  such  similarities. 

In  juxtaposition  with  the  formulation  and  implementation  of  the  simulation  code, 
we  have  devoted  significant  efforts  to  the  development  of  relevant  methodologies  for 
analysis  of  the  “data”  generated  by  the  molecular  dynamics  simulations.  As  an  example, 
one  important  program  that  we  have  developed  allows  us  to  go  through  the  data  of  the 
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entire  simulation  and  set  virtual  “detection  planes”  (see  Fig.3)  normal  to  the  z-axis 
(electric  field  direction)  at  specified  fixed  distances  above  and  below  the  parent  drop,  and 
record  the  statistics  of  all  clusters  that  pass  across  the  planes.  Detailed  information  on 
these  clusters  is  written  to  data  files  fi-om  which  a  variety  of  statistics  can  be  quickly 
extracted.  One  also  has  the  information  saved  so  that  any  of  the  emitted  clusters  can  be 
taken  and  allowed  to  continue  evolving  individually  in  a  separate  simulation. 

To  illustrate  some  of  the  more  easily  deducible  statistics,  we  show  in  Fig.6 
information  derived  from  the  simulation  discussed  above,  i.e.  a  formamide  droplet, 
containing  29  %  (by  weight)  of  dissolved  Nal,  placed  in  a  uniform  1.5  V/nm  external 
field.  The  detection  planes  were  placed  at  +200  nm  along  the  z-axis  (see  Fig.3).  The 
figure  shows  the  abundances  of  overall  flow  statistics  derived  from  the  clusters  passing 
these  detector  planes.  In  a  simulation  one  is  able  to  look  in  considerable  detail  at  what 
contributes  to  the  flow  of  material  emanating  from  the  ends  of  the  elongated  droplet. 
Fig.7(a)  depicts  the  abundances  of  ion-containing  clusters  of  different  charge  and  shows 
that  most  of  the  observed  clusters  have  charge  ±1.  A  decomposition  of  this  information  is 
seen  in  Fig.7(b)  that  shows  that  in  these  abundant  clusters  (charge  +1)  the  largest 
contributions  are  from  single  ion  clusters  (Na^  and  I  with  adsorbed  formamide 
molecules). 

More  detailed  information  is  displayed  in  Fig.8(a)  showing  die  number  of  singly 
charged  clusters  observed  (i.e.  positive  Na^,  Na2’^  T,  etc.  and  negative  T,  Na^  l2~,  etc.) 
without  regard  to  the  number  of  accompanying  adsorbed  formamide  molecules.  While 
multi-ion  clusters  do  contribute,  the  primary  contribution  is  from  single-ion  clusters  Na"^ 
and  r  and  one  sees  that  the  ratio  of  Na^  to  F  clusters  is  ~  2:1.  Fig.8(b)  provides 
information  related  to  the  size  and  mass  distribution  of  these  single-ion  clusters  by 
showing  the  relative  number  of  single-ion  clusters  having  a  specific  number  of  adsorbed 
formamide  molecules.  The  Na"^  clusters  have  predominantly  4-5  formamides  while  the  I" 
clusters  exhibit  a  much  broader  distribution.  It  is  interesting  to  note  that  experimental 
estimates  (utilizing  colloid  thrusters  and  formamide/Nal  mixtures)  of  the  average  number 
of  adsorbed  formamides  lay  within  in  this  range  of  4-5.  This  is  also  consistent  with 
spectroscopic  analysis  (ref  13)  of  Na(FA)n^  clusters  (n  is  the  number  of  formamide 
molecules  HCONH2,  denoted  here  by  FA).  This  experimental  work  employs  a  formamide 
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fluid  with  the  same  concentration  of  Nal  as  used  in  our.  simulations  and  instead  of  a 
droplet  in  an  external  field,  it  used  a  colloid  thruster  (at  the  mm  length  scale)  to  produce 
the  charged  clusters.  The  spectroscopic  analysis  indicates  that  for  n  >  4  the  abundance  of 
clusters  falls  off  dramatically.  This  is  of  interest  because  if  one  estimates  the  energy 
necessary  to  remove  one  formamide  molecule  from  Na(FA)n^  (see  Fig.9,  the  energies 
displayed  were  obtained  fi'om  simulations  we  performed  of  Na(FA)n'^  at  300K  for  n=l-8) 
there  is  no  abrupt  change  as  n  exceeds  four.  This  is  also  consistent  with  other  theoretical 
work  (ref  13)  where  the  energies  were  derived  from  optimized  Na(FA)n^  structures  and 
quantum  chemical  calculations.  Therefore,  it  is  rather  interesting  and  promising  that 
results  fi'om  our  simulations,  involving  drops  of  FA/Nal  mixtures  at  ~10^  nm  length 
scales  and  V/nm  applied  electric  field  strength,  already  begin  to  display  similar 
phenomena  and  statistics  as  found  in  mm  scale  colloid  thrusters  where  the  field  strengths 
due  to  the  applied  nozzle/extractor  potentials  can  be  considerably  smaller. 

We  have  found  that  if  the  detection  planes  are  moved  closer  to  the  parent  cluster, 
e.g.  ±100  nm,  the  formamide  distribution  (such  as  the  one  shown  in  Fig.  8(b)  for 
detection  planes  that  are  placed  at  ±200  nm)  shifts  very  slightly  to  higher  numbers.  This 
leads  us  to  conclude  that  the  effects  of  molecular  evaporation  may  be  observed  even 
relatively  close  to  the  emission  source.  This  touches  on  an  important  point  regarding 
these  types  of  simulations  and  dieir  relation  to  colloid  thruster  research.  In  laboratory 
experiments  the  charged  clusters  are  studied  far  ‘downstream’  firom  the  emitting  source 
and  the  statistics  of  the  clusters  that  are  observed  may  be  quite  different  than  just  after  the 
clusters  leave  the  emitting  jet  since  a  very  large  amount  of  evaporation,  secondary 
emissions  and  fragmentation  may  take  place.  The  mechanism  of  any  thrust  produced  by 
the  emitted  charged  species  depends  on  its  charge  and  mass  starting  at  the  instant  of 
leaving  the  source,  whether  it  is  a  Taylor-cone  or  parent  droplet.  Since  these  simulations 
allow  us  to  study  in  great  detail  the  natvire  of  the  charged  species  as  they  are  actually 
being  produced,  we  are  in  a  unique  position  to  provide  useful  insights  on  propellant 
properties  and  their  role  in  thrust  generation. 

The  work  discussed  above  serves  as  a  prelude  to  simulations  of  a  colloid  duuster 
where  one  wishes  to  understand  in  some  detail  the  factors  contributing  to  the  thrust.  The 
center-of-mass  of  the  cluster  in  the  simulation  discussed  above  is  stationary  so  that  the 
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simulation  may  be  regarded  as  mimicking  two  colloid  thrusters  back-to-back  with 
opposite  nozzle/extractor  bias  voltages.  The  thrust  displayed  in  Fig.6  is  the  average  of  the 
two  thrusts  (their  absolute  values)  due  to  charged  clusters  being  pulled  from  the  droplet 
and  accelerated  along  the  ±  z-axis,  parallel  to  the  imposed  field.  Referring  to  Fig.  10  we 
will  discuss  one  example  of  the  analysis  tools  that  we  have  explored  to  better  understand 
contributions  to  thrust  from  different  types  of  charged  clusters  (the  thrust  calculations 
here  are  due  to  the  positively  charged  clusters  only).  One  way  to  define  a  per-cluster 
thrust  contribution,  T(i),  is  to  define  a  local  mass  flow  rate,  associated  with  a  cluster,  e.g. 
file  i-th  cluster,  to  use  in  T(i)=dm/dt  x  u(i).  The  average  of  T(i)  over  all  of  the  clusters 
gives  an  average  thrust.  Another  way  to  estimate  the  overall  thrust  consists  of  evaluation 
of  the  quantity:  (total  average  mass  flow  rate)  x  (avg.  flow  speed).  Although  the  flow  is 
highly  non-uniform  and  composed  of  many  cluster  types,  both  approaches  are  found  to 
yield  the  same  order  of  magnitude. 

Fig.  10  displays  the  per-cluster  contribution  as  a  function  of  the  total  number  of 
ions  (positively  and  negatively  charged)  in  the  cluster.  Note  that  clusters  with  a  large 
number  of  ions  and  large  mass  have  high  individual  thrust  but  there  are  relatively  few  of 
them  so  they  don’t  contribute  significantly  to  the  average  thrust,  while  clusters  having 
few  ions  and  contributing  less  (individually)  to  the  total  thrust  are  heavily  weighted  by 
their  large  numbers,  and  thus  they  comprise  the  largest  tigaist  component.  The  idea  is  that 
individual  clustCT  thrusts  as  well  as  energies,  etc.  can  be  correlated  to  a  clusters  mass, 
charge,  number  of  ions,  etc. 

Phase  Two:  Complex  Geometries  and  Electric  Fields 

While  some  of  the  issues  pertaining  to  the  properties  of  liquid  droplets  in  electric  fields 
will  be  pursued  fiirther,  this  work  allowed  us  to  proceed  towards  our  primary  goal  -i.e., 
simulations  of  electrified  nanojets  and  colloid  thrusters.  In  the  second  phase  of  the 
project  the  computer  code  constructed  for  the  purpose  of  simulations  of  droplets  formed 
the  basis  for  the  development  of  a  code  that  can  model  the  complex  geometry  of  a 
capillary  needle  and  counter  electrode  (extractor  plate)  of  a  colloid  thruster.  In  particular. 
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the  electric  field  produced  by  an  applied  potential  difference  between  the  needle  and  plate 
is  quite  complex  and  this  must  be  integrated  into  the  parallel  fast  multipole  method 
(FMM)  code  that  was  described  earlier.  There  are  also  several  additional  components  that 
are  required.  For  example:  (i)  there  has  to  be  a  means  of  replenishing  the  fluid  atoms  lost 
fi-om  the  capillary  as  they  exit  and  are  accelerated  by  the  external  electric  field;  (ii)  one 
has  to  be  able  to  feed  a  fluid  (in  this  case  a  dielectric  fluid  with  dissolved  salts)  throu^ 
the  capillary  needle  either  at  a  prescribed  velocity  or  by  means  of  a  constant  applied 
pressure;  and  (iii)  an  algorithm  must  be  provided  for  neutralizing  the  charged  ionic 
species  that  remain  within  the  capillary  as  charge  is  transported  away.  In  addition,  we 
have  made  an  effort  to  assure  that  the  implanentation  of  these  elements  would  reflect,  as 
closely  as  possible,  the  associated  physical  mechanisms  present  in  an  actual  thruster. 

The  simulation  is  performed  using  the  modified  FMM  code  described  earlier, 
where  a  basic  cubic  root-cell  is  duplicated  along  one  axis  (the  z-axis)  to  create  a  region  of 
space  that  is  sufficiently  large  so  that  the  design  of  the  material  elements  and  tiie  required 
fields  will  fit  within  it;  this  is  coupled  of  course  to  the  requirement  that  the  long-range 
forces  of  the  charged  species  will  be  computed  accurately.  The  focus  is  on  colloid 
thrusters  which  typically  involve  a  metallic  capillary  tube,  or  nozzle,  possibly  with  a 
conical  tip  to  focus  the  fields  at  it’s  tip,  and  a  extractor  plate,  with  a  large  hole,  placed  a 
certain  distance  away  fi-om  the  nozzle.  The  extractor  plate  is  maintained  at  a  chosen 
potential  relative  to  the  nozzle  (see  Fig.l  1). 

Reservoir.  To  address  the  issue  of  a  reservoir  that  will  continually  supply  the 
nozzle  with  dynamic  fluid  atoms,  we  apply  techniques  that  have  been  used  in  computer 
simulations  of  the  formation,  stability  and  breakup  of  nanojets  (ref  14),  modified  to 
perform  with  an  ionic  fluid  (employing  the  FMM  structure  of  our  code),  and  designed  to 
allow  an  effectively  unlimited  amount  of  reservoir  material.  A  bulk  fluid  of  the  desired 
material  (formamide  with  a  prescribed  concentration  of  Nal)  is  prepared  at  ttie  desired 
working  temperature  (usually  ~  300K)  and  pressure.  The  first  step  is  to  take  a  cylindrical 
core  of  fluid  (fi:om  the  prepared  bulk  fluid)  that  will  be  placed  and  positioned  inside  the 
capillary  tube.  To  do  this,  copies  of  the  3-D  periodic  calculation  cell,  used  in  the  bulk 
fluid  calculation,  are  placed  together  so  that  the  width  of  the  resulting  block  of  fluid  is 
larger  than  file  desired  nozzle  inner  diameter,  and  the  length  fills  the  nozzle.  The  fluid 
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atoms  outside  the  desired  radius  are  discarded  (under  the  constraint  of  overall  charge 
neutrality).  In  this  way  a  relatively  small  bulk  sample  can  be  equilibrated  and  used  to  fill 
a  nozzle.  This  fluid  core  still  has  a  number  of  periodically  replicated  cells  along  the  z- 
axis.  A  boundary  is  defined,  acting  effectively  as  a  “piston”,  above  which  fluid  atoms  are 
treated  dynamically  and  below  which  the  reservoir  fluid  atoms  are  effectively  frozen  and 
the  forces  between  these  atoms  need  not  be  computed  thereby  saving  a  considerably 
amount  of  computational  time  (although  forces  between  these  atoms  and  the  dynamic 
atoms  are  computed).  The  piston’s  z-axis  position  is  exactly  one  leaf  cell  width  (the 
smallest  cell  size  of  the  FMM  tree  structure,  see  earlier  discussion)  above  the  bottom  of 
the  FMM  root-cell  -  this  is  purposefully  designed  so  that  the  dynamic  atoms  plus  one 
layer  of  cells  containing  static  reservoir  atoms  are  within  the  root-cell  (where  forces  can 
be  computed),  with  the  rest  of  the  reservoir  atoms  being  located  outside  the  root-cell. 
This  allows  all  of  the  dynamic  atoms  to  interact  with  a  portion  of  the  static  atoms. 

The  static  reservoir  atoms  are  moved  (or  pushed)  forward  into  the  lower  end  of 
the  nozzle  (the  method  of  moving  is  discussed  below).  All  of  the  atoms  in  the  system  fall 
into  the  categories  of  dynamic  or  static.  As  the  static  reservoir  atoms  cross  the  piston  z- 
position  they  are  treated  fully  dynamically.  In  the  case  of  molecules,  when  a  molecule’s 
center-of-mass  (cm)  crosses  the  piston  static/dynamic  boundary  the  molecule’s  atoms  are 
declared  d3mamic.  Once  an  atom/molecule  is  declared  as  d)mamic,  it  remains  in  this 
category  for  the  duration  of  the  simulation.  For  a  prescribed  distance  above  the  piston 
(we  use  3  nm),  the  newly  dynamic  fluid  atoms  are  temperature  controlled  via  stochastic 
thermalization.  After  this  point  the  fluid  atoms  that  are  within  a  chosen  distance  of  the 
nozzle  walls  are  temperature  controlled. 

To  discuss  the  method  by  which  the  reservoir  material  is  replenished,  recall  that 
the  cylindrical  nozzle/reservoir  material  is  comprised  of  several  rqjlieations,  n,  of  a 
periodic  system  with,  e.g.,  a  cell  of  width  Wz  along  the  z-axis.  Whenever  an  atom,  or  a 
molecule  (i.e.  the  molecule’s  center-of-mass,  cm),  crosses  flie  piston  boundary  and 
becomes  dynamic,  an  exact  copy  of  it  is  placed  at  the  rear  of  the  advancing  “frozen” 
reservoir  by  subtracting  nWz  from  the  z-coordinate  of  the  particle,  and  adjusting  the 
result  for  the  past  advaneements  of  the  reservoir’s  cm.  Using  this  method,  one  can 
construct  the  nozzle/reservoir  fluid  from  a  few  replications  of  a  relatively  small  unit  cell 
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and  reuse  the  reservoir  material  over  and  over.  As  mentioned,  just  beyond  the  piston,  the 
fluid  molecules  are  stochastically  thermalized,  erasing  any  effects  of  periodicity. 

To  move,  or  push,  the  reservoir  material  forward  through  the  nozzle,  we  use  one 
of  two  relatively  simple  methods.  In  the  first,  which  we  are  currently  using,  the  reservoir 
material  is  moved  forward  at  a  constant  velocity.  The  second  method  corresponds  to  a 
constant  pressure  being  applied  to  the  reservoir  fluid.  The  forces  on  the  reservoir  atoms 
with  flie  dynamic  atoms  across  the  piston  boundary  are  computed  (divided  by  the  piston’s 
area,  they  give  the  time  varying  pressure  on  the  piston  and  reservoir  material).  The 
desired  constant  externally  applied  pressure  is  multiplied  by  the  piston  area  to  produce  its 
associated  force.  The  difference  between  this  force  and  the  calculated  force,  along  with 
the  current  cm  velocity  of  the  reservoir  core,  is  then  used  in  a  simple  integration 
algorithm  to  advance  all  of  the  reservoir  material  by  a  common  amount  at  each  time  step. 

Treating  excess  charges.  In  a  colloid  thruster,  charged  droplets  are  created  and 
accelerated  through  the  extractor  plate  resulting  in  charge  imbalances  and  buildup  in  the 
thruster  unit.  So  if  the  working  fluid  is  formamide/Nal,  then  a  negatively  biased  extractor 
plate  will  extract  positively  charged  material  (excess  of  Na"^)  leaving  an  excess  of  T  ions 
at  the  nozzle.  It  is  currently  believed  that  these  ions  lose  their  charge  to  tiie  nozzle  and 
form  neutral  molecules  which  are  then  lost  through  several  processes,  including:  (i) 
evaporation,  (ii)  convective  flow  into  the  charged  material  leaving  the  system,  (iii) 
diffusion  either  on  the  exterior  the  nozzle  surface  or  back  into  the  nozzles  reservoir 
material.  Namely,  the  general  picture  is  that  the  excess  charged  ionic  species  left  behind 
loses  it’s  charge  somewhere  near  the  metallic  nozzle  and  flie  resulting  neutral  species  is 
effectively  removed,  without  effecting  the  operation  of  the  thruster.  To  model  this  mode 
of  operation,  and  maintain  a  charge  balance  in  the  simulation  of  a  colloid  thruster,  we 
have  created  the  following  algorithm.  First  a  cylindrical  test  section  is  specified,  having  a 
narrow  width  along  the  z-axis  that  can  be  varied  and  having  a  radius  that  exceeds  the 
expected  base  radius  of  the  fluid  before  it  tapers  down  due  to  the  action  of  the  externally 
applied  field.  The  net  charge  from  the  atoms  within  the  test  volume  is  computed  and 
when  the  charge  becomes  negative  (for  the  case  of  excess  T  ions)  the  outermost  T  ion 
near  the  periphery  of  the  test  volume  is  removed  from  the  system.  Note  that  this  does  not 
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mean  that  the  tapering  form  of  the  exit  fluid,  e.g.  a  Taylor  cone,  need  be  neutral.  In  fact, 
in  a  Taylor  cone  the  surface  has  a  charge  density,  and  it  is  the  interplay  between  capillary 
forces  and  the  surface  forces  due  to  the  external  field  that  permits  a  steady  state  conical 
shape.  Once  a  steady  state  is  reached  in  the  simulation,  one  expects  that  any  excess 
charge  beyond  the  surface  charge  will  be  neutralized  near  the  nozzle  and  the  above 
algorithm  assures  this  by  maintaining  neutrality  in  a  zone  near  the  nozzle  exit. 

Boundaries.  Before  discussing  the  implementation  of  the  external  field  we  give 
a  few  details  about  the  rules  concerning  the  motion  of  atoms  and  the  system  boundaries. 
Boundaries  of  the  system  include  the  boundaries  of  the  FMM  root-cell  (atoms  that  cross 
these  are  simply  omitted  fi-om  the  system)  as  well  as  the  metallic  nozzle/plate  material. 
The  nozzle  is  created  out  of  a  solid  block  of  a  chosen  metal  by  discarding  the  atoms 
outside  the  boundaries  that  define  the  nozzle.  The  extractor  plate  is  created  in  a  similar 
manner  by  using  a  defining  boundary  that  is  relatively  thin,  has  a  large  circular  hole,  and 
is  centered  along  the  z-axis  at  the  desired  distance  fi-om  the  nozzle  exit.  The  diameter  of 
the  hole  is  close  to,  and  is  smaller,  than  the  width  of  a  FMM  root-cell,  and  the  outer  x-y 
boundaries  of  the  extractor  plate  are  outside  of  the  FMM  root-cell  (The  geometry  and 
positioning  of  the  nozzle  and  plate-electrode  are  further  discussed  in  latter  sections). 

Particles  can  clearly  not  be  allowed  inside  the  nozzle/plate  material  and  this  is 
realized  though  the  repulsive  component  of  the  particle-wall  interactions.  Namely,  a 
primary  function  of  the  particle-wall  forces  is  to  keep  atoms  jfrom  crossing  wall- 
boundaries,  and  these  forces  can  vary  fi-om  purely  repulsive  forces,  modeling  a  non¬ 
wetting  fluid,  to  a  balance  of  attractive  and  repulsive  forces  (approximating  known  or 
desired  interactions).  These  particle  interactions  with  the  walls  (i.e.  nozzle  or  extractor 
plate  material)  can  be  implemented  in  a  number  of  ways.  One  method  is  to  simply  use 
appropriate  atomic  interaction  potentials  between  all  of  the  fluid  material  and  all  of  the 
nozzle/plate  atoms.  This  would  necessitate  a  substantial  increase  in  the  required 
computational  load.  The  basic  operation  of  the  simulated  thruster  is  expected  to  be 
insensitive  to  certain  details  of  these  interactions  and  one  can  explore  simpler 
representations  of  the  particle-wall  interactions  that  can  be  evaluated  in  a  highly  efficient 
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manner.  With  these  considerations  in  mind,  we  describe  below  Ihe  manner  in  which  we 
have  implemented  atom-wall  interactions. 

When  charged  (or  neutral)  material  encoxmters  the  extractor  plate  outside  of  its 
hole  radius,  it  may  be  expected  to  be  neutralized  (or  remain  neutral)  and  eventually 
migrate  outside  of  the  operational  area  of  the  thruster.  As  a  first  approximation  of  this 
circumstance,  when  any  atom  or  molecule  crosses  the  boundary  defining  the  extractor 
plate,  they  are  removed  from  the  system.  With  regard  to  atom-nozzle  interactions  we 
note  that  tiie  most  important  elements  of  the  thruster  operation  are  capillary  forces  and 
the  interactions  of  the  external  field  with  the  fluid  atoms  that  have  left  the  immediate 
vicinity  of  the  nozzle  exit.  Thus,  for  example,  if  the  assumed  atom/nozzle  forces  lead  to 
wetting  of  the  exterior  nozzle  tip  then  the  base  of  a  Taylor  cone  will  be  anchored  at  the 
nozzle  tip  (as  it  should),  and  this  wetting  action  can  be  achieved  throu^  forces  that  are  a 
simple  and  efficient  representation  of  more  complex  particle- wall  forces.  We  have  just 
entered  in  our  simulations  the  stage  of  actual  thruster  operation  and  in  line  with  the 
preceding  discussion  we  initially  take  the  nozzle  walls  to  be  smooth  in  the  sense  that 
interaction  potentials  are  taken  to  depend  only  on  tiie  distance  of  the  fluid  atoms  from  the 
wall.  For  atoms  close  to  a  wall,  we  model  the  atom- wall  interaction  as  a  3-9  Lennard- 
Jones  (LJ)  potential  (that  describes  the  interaction  of  an  atom  interacting  with  atoms 
comprising  a  half-space  solid  through  a  6-12  LJ  potential);  the  parameters  of  the 
interaction  potential  are  chosen  to  be  close  to  those  of  hydrocarbon-gold  interactions  (ref 
15).  However,  as  discussed  above,  the  actual  values  of  the  potential  parameters  are  not 
expected  to  be  of  great  importance,  and  they  may  be  varied  to  achieve  the  appropriate 
wetting  behavior. 

Electric  field.  At  this  point,  we  will  discuss  one  of  the  most  important  elements 
in  the  development  of  our  simulation  code,  namely  the  addition  of  the  electric  field  that 
depends  on  the  nozzle  and  plate  geometry  and  the  potential  difference  between  them. 
Once  the  nozzle/plate  geometry,  size,  and  placement  have  been  chosen,  we  input  them 
(along  with  the  required  potential  difference)  into  a  finite-element-method  (FEM) 
program  (see  ref.  16)  used  to  calculate  the  potential,  0(z,r),  and  electric  field,  E(z,r),  on  a 
cylindrically  symmetric  2-D  grid  that  is  large  enough  to  encompass  the  entire 
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calculational  cell  (the  root-cell).  For  example,  in  Fig.  12  we  show,  with  axial  symmetry, 
the  FEM  boundaries  (thick  lines)  that  are  positioned  to  lie  near  the  surfaces  of  the  nozzle 
and  extractor  plate.  A  higher  density  of  nodes  is  used  near  the  nozzle  tip  to  obtain  better 
resolution  of  the  electric  field. 

The  set  of  field  values  calculated  at  the  nodes  of  the  FEM  mesh  are  mapped  onto 
a  cylindrical  grid,  and  they  are  used  as  input  to  the  colloid  thruster  MD  simulation  code. 
At  any  given  time,  the  charged  ions,  or  the  molecular  components  having  partial  charges, 
have  cylindrical  coordinates  lying  within  one  of  the  grid  cells  on  whose  voices  the  field 
values  are  known.  Using  the  particle  coordinates  relative  to  the  grid  cell  vertices  one  can 
then  use  a  simple  2-D  interpolation  routine  to  compute,  in  cylindrical  coordinates,  the 
potentials  0(z,r)  and  fields  E(z,r)  at  the  location  of  each  of  the  charged  particles.  After 
resolving  the  electric  field  into  its  three  components  relative  to  the  orthogonal  Cartesian 
axes,  the  force  on  an  atom,  with  charge  q,  is  computed  as  F(r)  =  qE(r)  and  potential 
energy  U(r)  =  q4>(r)  where  r  =  (x,y,z). 

Further  details  about  the  colloid  thruster  simulations,  as  well  as  observations 
concerning  colloid  thruster  or  electrospray  operation  at  smaller  length  scales,  can  be 
discussed  by  reference  to  Fig.13,  where  we  display  a  cross  section  of  the  FMM  root-cell 
with  areas  that  bound  the  metallic  nozzle/extractor  material  shown  as  gray  solid  areas.  It 
is  within  these  areas  tiiat  the  FEM  boundaries  (Fig.  12)  are  prescribed.  The  variation  of 
the  electric  field  strength  is  shown  through  the  use  of  color  gradations  and  the  arrows 
point  in  the  field  direction.  Equipotential  lines  are  also  shown.  The  extractor  plate/nozzle 
potential  difference  is  —50  V  and  because  of  the  ~  10^  nm  length  scales,  the  field 
strengths  reach  values  that  are  of  the  order  of  V/nm  (see  the  Future  Plans  section  for  a 
discussion  of  these  field  strengths  as  compared  to  those  encountered  in  devices  defined 
on  the  mm  length-scale). 

Testing  of  the  nano-thruster  simulation  code.  Our  present  simulations,  using 
die  geometry  and  field  values  shown  in  Fig.13,  have  involved  non-wetting  walls.  A 
sample  of  our  results  is  shown  in  Figs.  14-16.  We  observe  fi-om  Fig.l4  (and  a  closer  view 
provided  in  Fig.  15),  that  the  fluid  cone  emanating  from  the  nozzle  extends  down  to  the 
piston;  in  this  simulation  the  fluid  is  pulled  by  the  applied  electric  field,  as  well  as  pushed 
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through  the  nozzle  at  a  velocity  of  6n3/s.  The  speed  of  flow  of  the  fluid  jet  in  a  3  nm-wide 
region  near  the  exit  is  50  m/s.  The  flow  accelerates  towards  the  extraction  plate,  reaching 
velocities  of  3000-5000m/s  at  and  just  past  the  plate  (with  the  lower  value  corresponding 
to  a  formamide-solvated  sodium  iodide  molecule,  and  the  higher  one  corresponding  to  a 
solvated  single  sodium  ion). 

Other  views  of  the  simulated  system  are  shown  in  Fig.  16.  This  describes  a  nozzle 
whose  wall  material  is  unwetted  by  the  fluid  above  a  certain  point,  while  below  this  point 
wetting  occurs.  In  this  study  the  two  regions  are  of  the  same  radius;  however,  we  have  the 
freedom  to  modify  this  configuration  to  one  where  the  radius  of  the  upper  exit  region  is 
larger  than  that  of  the  region  to  the  left  of  the  piston  (see  Fig.  15). 

The  material  streaming  from  the  leading  edge  of  the  exiting  fluid  consists  of  both 
evaporating  neutral  formamide  molecules  as  well  as  positively  charged  clusters  that  are 
accelerated  by  the  external  fields  through  the  extractor  plate.  Because  of  the  high 
acceleration  of  detaching  charged  clusters,  one  generally  observes  only  a  few  clusters 
between  the  end  of  the  fluid  cone  and  extractor  plate.  Once  past  the  extractor  plate,  the 
charged  clusters  are  in  a  much  weaker  field  and  they  move  essentially  at  a  constant 
velocity.  Note  flie  similarities  between  the  flow  of  the  detaching  charged  clusters  in 
Figs.  14-16  for  a  colloid  thruster  and  in  Figs.4-5  for  a  droplet,  in  line  with  the  expectation 
that  both  share  a  common  physical  description.  The  analysis  programs  used  to  calculate 
file  statistics  of  emitted  clusters  from  droplets  in  a  constant  electric  field  are  being 
modified  for  use  in  the  colloid  thruster  simulation.  For  example,  a  virtual  ‘detection 
plane’  (Fig.3)  may  be  placed  at  any  desired  distance  from  the  nozzle  exit  and  the  statistics 
of  the  material  crossing  this  plane  collected  and  analyzed  (the  region  of  primary  interest 
is  just  past  the  extractor  plate). 

Summary.  To  summarize  this  section,  we  remark  that  we  have  arrived  at  a  stage 
where  we  have  constructed,  implemented  and  tested  an  efficient  simulation  code  for 
investigations  of  the  properties  and  mechanisms  underlying  the  physical  characteristics 
and  operational  principles  of  a  colloid  thruster.  While  some  work  still  needs  to  be  done  to 
bring  the  simulation  model  closer  to  its  experimental  counterpart,  the  nontrivial  task  of 
constructing  the  most  essential  compute  algorithms  is  completed. 


19 


Acknowledgements:  We  thank  Rainer  Dressier  for  helpftil  discussions  and  for 
sharing  his  work  with  us. 


References: 

1)  M.  P.  Allen  and  D.  J.  Tildesley,  Computer  Simulation  of  Liquids  (Clarendon,  Oxford, 
1987) 

2)  M.  Svanberg,  Mol.  Phys.  92, 1085  (1997). 

3)  E.  D.  Stevens,  Acta  Cryst.  B34,  544  (1 978). 

4)  W.  D  .Cornell,  et  al.,  J.  Am.  Chem.  Soc.  117, 5179  (1995). 

5)  www.amber.ucsf.edu/amber/ff94/parm.dat 

6)  E.  Sigfridsson  and  U.  Ryde,  J.  Comp.  Chem.  19,  377  (1998). 

7)  L.  Perena,  et  al.,  J.  Chem.  Phys.  102, 450  (1995). 

8)  L.  F.  Greengard,  The  Rapid  Evaluation  of  Potential  Fields  in  Particle 
Systems  (The  MIT  Press,  1988). 

9)  The  physics  is  not  expected  to  be  significantly  different  than  for  T=300K. 
Nevertheless,  we  work  at  a  slightly  higher  temperature  31  OK  to  account  for  the  10  K 
higher  melting  point  of  the  simulated  formamide  so  that  our  results  correspond  more 
closely  to  300K  experimental  results. 

10)  M.  Gamero-Castano  and  V.  Hruby,  J.  of  Propulsion  and  Power  17, 977  (2001). 

1 1)  M.  Gamero-Castano  and  J.  Fernandez  de  la  Mora,  J.  Chem.  Phys.  113,  815  (2000). 

12)  H.  A.  Stone,  J.  R.  Lister,  and  M.  P.  Brenner,  Proc.  R.  Soc.  London  A 
455,329(1999). 

13)  Rainer  Dressier,  private  commxmication.  Also  a  presentation  by  Yu-Hui  Chiu  of  his 
group  at  the  Colloid  Thruster  /  Nano  Electrojet  Workshop,  Cambridge,  MA,  Oct.  3- 
4,2000. 

14)  M.  Moseler  and  U.  Landman,  Science  289,  1165  (2000). 

15)  T.  K.  Xia,  J.  Ouyang,  M.  W.  Ribarsky  and  U.  Landman,  Phys.  Rev. 
Lett.  289, 69  (1992). 


20 


16)  We  use  a  readily  available  FEM  code  that  may  be  freely  downloaded  at 
www.quickfield.com.  This  code  is  simple  to  use  so  that  changes  in  the  nozzle/extractor 
plate  design  and  desired  boundary  conditions  are  relatively  easy  to  implement. 

1 7)  J.  Fernandez  de  la  Mora,  private  communication. 

18)  J.  Fernandez  de  la  Mora  and  I.  G.  Loscertales,  J.  Fluid  Mech.  260, 155  (1994). 

19)  I.  G.  Loscertales,  and  J.  Fernandez  de  la  Mora,  J.  Chem.  Phys.  103,  5041  (1995). 


21 


FIGURES 


Effect  of  an  electric  field  on  a  formamide  droplet 
molecular  dynamics  simulation 
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a  =50  mn 


Fast  Multipole  Method 


•  Binary  Tree  Structure  Division  of  Space 
^  Atoms  in  adjacent  cells  interact  directly 

•  Compute  multipole  moments  of  all  cells  at  all  levels 
from  the  atoms  within  the  cells 

•  Atoms  interact  only  with  the  multipoles  of  more  distant  cells 


Fig  2 
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153.6  nm 


System  Geometry 

‘detector  planes’  £ 


460.8  mn 


•  7150  Formamide  molecules,  810  Na^,  8101 

(detector  planes  normally  set  at  ±  200  nm)  Fig  3 
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MD  simulations  of  electric-field-induced  fragmentation  of 
aformamide  /Nal  droplet  (29%  Nal  by  weight) 
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Cluster  Flow  Statistics 


averaged  over  +  and  -  flow 

_ E=1.5V/nm _ 

Total  flow  rate  ~  4x10*^  nl/s 


Charged  flow  rate  ~  2  x  10"^  nl/s 

Current  ~  30  nA 
Thrust~  1.6  nN 


Fig  6 
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Cluster  Charge  Statistics  E=1.5  V/nm 

(for  clusters  having  an  ionic  component) 


Cluster  charge 


Fig  7  (a) 
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Cluster  Charge  Statistics  E=  1 . 5  V/nm 

(for  clusters  having  an  ionic  component) 
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Singly  Charged  Cluster  Statistics  E=1 .5  V/nm 
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Exploring  ways  of  understanding  results:  thrust  contributions  Tj 
tj  =  time  that  clusterj  passed  detector  plane 

rnj.Uj 


Ti  =  (m|/Tj)xUi 

T  =  avg(T|) 

~1 .6  nN 


Note: 
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Extractor  Plate 


Positively  charged  clusters 
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Colloid  Thruster  MD  Simulation  Geometry  and  Fields 
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Colloid  Thruster  MD  Simulation 
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II.  Major  Accomplishments 

Formulation,  development,  implementation  and  testing  of  efficient  and  flexible 
computational  methodologies,  algorithms  and  parallel  computer  codes  for  simulations 
and  analysis  of : 

(a)  Equilibrium  and  non-equilibrium  properties  of  neutral  or  charged  fluid  droplets 
(consisting  of  up  to  tens  of  thousands  of  atoms),  under  field-fi-ee  conditions  and  in  the 
presence  of  applied  electric  fields. 

(b)  Generation,  stability  and  breakup  of  fluid  jets,  driven  by  applied  electric  fields  - 
modeling  colloid  thrusters  of  nanoscale  dimensions. 

The  above  simulation  and  analysis  tools  have  been  applied  to  investigations  of 
geometrical  and  energetic  properties  of  droplets  and  nanojets,  made  of  a  salt  (Nal,  29% 
by  weight)  dissolved  in  formamide.  (see  section  (I)  of  the  report). 

III.  Interactions  and  Technology  Transitions 

In  October  2002,  Uzi  Landman  (PI)  and  W.  David  Luedtke  presented  the  results  of  their 
studies  in  Boston  at  a  workshop  on  electrified  nanojets,  organized  by  Dr.  Rainer  A. 
Dressier  (Air  Force  Res.  Laboratory,  Space  Vehicle  Directorate,  Hanscom  AFB,  MA). 
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V.  Publications  and  Presentations 

Several  publications  pertaining  to  methodological  issues  and  description  of  the  physical 
results  of  our  simulations  of  dielectric  dipolar  droplets  (with  and  without  solvated  ionic 
salts)  in  strong  electric  fields,  are  in  the  process  of  preparation. 

Presentations: 

The  results  of  our  investigations  were  reported  at  the  March  2003  meeting  of  the  APS 
(American  Physical  Society)  in  Austin,  Texas,  and  the  most  recent  results  will  be 
presented  at  the  APS  in  Montreal  (March,  2004) 

Plenary  Speaker  (Landman),  Nano-7/Ecoss-21,  Malmo,  Sweden,  June,  2002. 

Speaker  (  Landman),  AFOSR  Workshop  on  Electrified  Nanojet  and  Colloidal  Thrusters, 
MIT,  Boston,  October,  2002.  Uzi  Landman  and  David  Luedtke  (a  Senior  Research 
Scientist  in  Landman’s  group  have  participated  in  the  workshop). 

Invited  Speaker  (Landman),  Conference  on  Cluster  and  Nanoscale  Science,  Nanjing 
University,  Nanjing,  China,  March,  2003. 

Keynote  Speaker  (Landman),  Conference  on  Challenges  in  Computational  Chemistry, 
Imperial  College,  London,  England,  April,  2003. 

Invited  Speaker  (Landman),  International  Union  of  Vacum  Societies  Conference  on 
NanoScience  and  Technology,  Eilat,  Israel,  May  2003. 
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